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Abstract 

In a tokamak pedestal, radial scale lengths can become comparable to the ion orbit width, 
invalidating conventional neoclassical calculations of flow and bootstrap current. In this work we 
illustrate a non-local approach that allows strong radial density variation while maintaining small 
departures from a Maxwellian distribution. Non-local effects alter the magnitude and poloidal 
variation of the flow and current. The approach is implemented in a new global 5f continuum 
code using the full linearized Fokker-Planck collision operator. Arbitrary collisionality and aspect 
ratio are allowed as long as the poloidal magnetic field is small compared to the total magnetic 
field. Strong radial electric fields, sufficient to electrostatically confine the ions, are also included. 
These effects may be important to consider in any comparison between experimental pedestal flow 
measurements and theory. 
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In the H-mode edge pedestal of a tokamak, strong density and temperature gradients 
drive large a neoclassical flow and bootstrap current. This flow and current affect stability 



of the region to ELMs and other modes. However, conventional 



neoclassical calculations 



are invalid in the pedestal since they rely on an expansion in the smallness of the 
poloidal ion gyroradius po to the perpendicular scale length of density and temperature r±. 
In the pedestal, this ratio pe/r± is not small. (We do not claim r± scales with p@, only 
that the lengths happen to be comparable in existing devices.) Physically, conventional 
neoclassical theory is based upon the smallness of the orbit width (~ pe for ions) relative 
to equilibrium profiles, yielding a local theory: flows and fluxes on one flux surface are 
determined by values and gradients of pressure p and temperature T and the electric field 
at that flux surface only. In the pedestal, however, equilibrium profiles can vary strongly on 
scale of the ion orbit width, requiring a global (nonlocal) calculation that does not rely on 
the conventional pe/f± expansion. 

In this work, we generalize neoclassical calculations both analytically and numerically 
to the case of a strong density pedestal (with density scale-length r n ~ pe) as long as the 
ion temperature scale length tt remains ^> pe, with a few other assumptions. We demon- 
strate how the neoclassical flow is altered, and the resulting poloidal flow variation will be 
important to consider for understanding experimental pedestal flow measurements. More 
generally, we emphasize that compared to the general r-r ~ Pe case, this "weak-T]' pedestal" 
is much more amenable to analysis: the distribution function remains nearly Maxwellian, 
permitting a Sf rather than full-/ approach and linearized treatment of collisions, and the 
E x .B-drift nonlinearity also becomes negligible. Any more ambitious effort to analyze a 
pedestal with ~ pg for finite aspect ratio will likely need to retain both these nonlineari- 
ties, necessitating complicated codes, which our results may be used to benchmark. We also 
present a new numerical continuum approach to computing these global neoclassical effects 
in the weak-TJ' limit, including the exact linearized Fokker-Planck collision operator. We 
exploit the success of local continuum neoclassical codes by making such a local code the 
inner step of an iteration loop for the global calculation. 

Several local neoclassical codes have been developed jsB], and other numerical efforts 
have computed nonlocal neoclassical effects in transport barriers using the particle-in-cell 



(PIC) approach 9Ml2j. Since PIC and continuum codes have differing treatments of col- 



lisions and boundary conditions and differing numerical resolution challenges, it is good 



practice to develop both approaches to ensure they yield the same physical results. Some 
neoclassical investigations have been made in global continuum codes [13I, Q], but these 
codes are ultimately designed for turbulence studies, and very different algori thms have 
been used than the one we describe. Some analytic results are available |l5l-ll 71]. but only 
in restricted limits of aspect ratio and collisionality, where simplified collision models are 
expected to be valid. 

Throughout our analysis we assume Bg <C B, where B = \B\is the magnetic field strength 
and Bg is the poloidal field, implying a scale separation between pg and the gyroradius p. 
Without this approximation, a gyrokinetic rather than drift-kinetic treatment would be 
necessary, including changes to the collision operator ; 

In conventional neoclassical theory, the ion distribution function is expanded as j\ = 
fm + fi where fui 3> fi, and /mi is a Maxwellian with constant density n\ and temperature 
T[ on each flux surface. The drift-kinetic equation is then solved for fi, with the result that 
/1 includes a term —(Iv\\/Q)dfMi/dijj. Here, / equals the major radius R times the toroidal 
field -Btor, ^ = ZeB/rriiC, Z is the ion charge in units of the proton charge e, m; is the ion 
mass, c is the speed of light, and 2irip is the poloidal flux. The derivative is carried out at 
fixed total unperturbed energy Wq = m- x v 2 /2 + Ze& , where $0 — is the flux-surface 
average of the electrostatic potential $. We may estimate dfui/dip ~ fMi/(RBgr±), so 
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fi ~ (pe/ r ±)fMi where pg = BvJ(Bgil) is the poloidal ion gyroradius, and V\ = ^2T i /m i 
is the ion thermal speed. In a pedestal, since pg/r± ~ 1, then fi ~ fuu so conventional 
neoclassical results are no longer valid. 



However, a more precise analysis reveals 15| a regime in which the near-Maxwellian as- 



sumption is still appropriate. Writing / M i = ^WWi/P^T^)]} 3 ^ 2 exp(— Wq/T-^)), where 
77(V>) = ni('0) exp (Ze$o(V , )/^ 1 i(V'))) the derivative (dfMi/dip)w that determines the magni- 
tude of fi is 

1 dr) ( 3\ 1 dTi] 



df M i 



dtp 

The magnitude of dfyii/dip is evidently determined by Tt and r v , the scale-lengths of T\ and 
rj, but not directly by r n , the scale-length of density. Observing r~ l = r" 1 — Ze^/Ti + 
Ze& / (Tir T ), }\/ fui m ay be small even when r n ~ p e as long as r T and r v are ^> pg. Such 
is the case when d^o/dip ~ T i (Zeni)~ l dn i /dip so the ions are electrostatically confined. 

We consider this "weak- T( pedestal" regime for the rest of the analysis: r n ~ pg but 
5^1 where 5 = pg/r T is the basic expansion parameter, and pg/r v ~ 5. (The electron 



temperature T e is free to vary on the pe scale.) This ordering, also considered in Ref. [15|, is 
useful in part because the collision operator may be linearized. Also, as we will show, the 
poloidal electric field decouples from the kinetic equation, so the equation becomes linear in 
fi. For Tt ~ Pe and/or r v ~ p#, the full bilinear collision operator must be used and a full-/ 
nonlinear kinetic equation must be solved, including the electric field nonlinearity. Notice 
r v <C Pe implies (Ze/Ti)d§o/dip ~ l/(RBgpg) and so Ze$o/T\ ~ 1. As a result, the term 
v e ■ V/i in the kinetic equation, neglected in conventional theory, becomes comparable in 
magnitude to the U||V||/i term. Thus, even though the weak-T/ ordering permits fi <C Jmu 
conventional neoclassical results still must be modified. As Bg <C B, the E x B drift Ve 
satisfies \v E \ Cw; so centrifugal effects may be neglected. 



We begin with the ion drift-kinetic equation 20] 



( W|| 6 + i; d )-(V/ i ) /X)W = C i {/i} + 5 (2) 

where the gradients hold fixed p = rriiV 2 / (2B) and W = m\v 2 /2 + Ze$ (now including $, 
not just $ )j Ci is the ion-ion collision operator linearized about fuu and S represents any 
sources/sinks. We take v& = (v\\/Q) V| w x (v\\b) (which includes ve-) 

Now change from W to Wo = W — Ze$i as an independent variable, where $i = $ — $o- 
We assume $i ~ <5<3>o and d<&\/dip ~ 5d§o/dijj, and we will show in a moment these orderings 
are self-consistent. Then defining g by 

fi = fm - (Ze^/TJfa - (I Vll /n)df Mi /d4> + g, (3) 

(j2J) may be written 

(v }] b + v d ) ■ iyg)^ Wo - dig} = d + S (4) 

where C\ = Ci{(Iv\\/Q)dfMi/dip} is the inhomogeneity, the independent variable is now 
Wo, and terms small in 5 have been dropped. The contribution from $! to v& ■ V9 is 
0(5) smaller than the $o contribution, and (v E • V^)/(u m • V?/>) ~ Ze^\jT\ ~ 5 where 
v m — v d — Ve is the magnetic drift, so we may approximate in (J3J) with the leading-order 
drift v d0 = v m + v E o where v E o = (c/B 2 )B x V$o- Then (T4]) is completely linear. To 
evaluate $i we may use the electron density n e + (e$i/T e )n e with quasineutrality to find 



e$i/T ; = (T-JT C + Z) 1 rz i 1 f d?v g. Hence, as g ~ Sfuu our assumed ordering for $! i 



is 



self-consistent. Using Ci{v\\f Mi} = 0, the only gradient surviving in C\ is dT/dtp. While the 



independence of g from dn-Jdip and d<&o/dip was known previously for the local case, the 
persistence of this property in the weak- 7]' pedestal case is noteworthy [l5| ■ 

One crucial difference between the local and global analyses is that the flow may vary 
over a flux surface in different ways. First consider the parallel ion flow niViii = J d 3 vv\\fi. 

cl {dpi d% B 2 dT\\ 

where ku = Ze (B 2 ) (clniB dTi/dip) -1 J d 3 vv\\g is dimensionless. We have exploited the 
aforementioned fact g oc dTi/dip. In the conventional ordering, fen is also the coefficient of the 
poloidal flow Vg: forming the appropriate linear combination of (JSJ) with the perpendicular 
diamagnetic and E x B flows, Vg = V ■ eg = k^cIBg (Ze (B 2 )) 1 dT\/dip where eg = 
(VC x VV>)/|VC x W| and Bg = B ■ eg. 

Applying J d 3 v = 2txty\ 2 J2a a I dW I dfJ>(B/v\\) to (@|, where a = sgn(u||), the resulting 
mass conservation equation (ignoring S) is 

9 f ^ t l \9 9 /" 3 gv m ■ Vji 
89 J dv ^ +u) B-d^j dv ^^ = ° 

where u = (vao ■ V6)fV\\6 w (cl / B)d<&o/dip is comparable in magnitude to In the local 
case, where v& ■ Vg is neglected in (j4]), only the first term in ((6]) (oc v\\) arises, implying 
J d 3 v vug oc B and dku/06 = 0. This is the origin of the well known conventional result that 
ku is constant on a flux surface. However, in the global case, the strong poloidal drift and 
pg-scale radial variation drive poloidal variation in ku. 

The total flow remains divergence-free in a fluid picture: V • T = where T = T\\b + TE + 
Tciia, r*] ] = / d 3 vv\\fi, Te = ti\Ve\ + / d 3 v fiVEo contains the first two orders of the E x B 
flux, v E i = (c/B 2 )B x V$i, r d ia = c(ZeB 2 )~ l B x V contains the first two orders of the 
diamagnetic flow, = Rl( I - bb) + p\\bb, pj_ = mj d 3 v fiv\/2, and p\\ = m J d 3 v fcv 2 . 
To prove V ■ T = from (jSJ), © an d Vnj ~ — (Zen i /T i )V < l ) o are applied, along with 
/ d 3 v fiv m = r dia + V x M + T f (true for any /■). Here M = bcp±/(ZeB), and we will 
neglect the 0(05) parallel flow correction Tf = (pu — p±)bb ■ V x b (which disappears when 
a more accurate v m is used.) As before, ni(ijj) = J d 3 v fyii includes only the leading-order 
density. We have needed to keep terms of two orders in both Te and T^a because the E x B 
and diamagnetic flows cancel to leading order in our ordering. And, though Fe ~ ti-iVeo 
and fa pi I , the radial derivative in V • T means the next-order corrections to these 
terms must be retained to accurately compute V • T. The poloidal fluid velocity is defined 



by Vg = r • eg/rii. It can be shown that to leading order in 8 



mB 



Cld%\ 18 r 3 v{ 

dv ^^B^) 9+ mJ dv f g . 



(7) 



In the local case, the v\\ term dominates, so Vg oc h\. In the global case, Vg remains 
proportional to dTjdip, but Vg is no longer oc k\\. A normalized poloidal flow may be 
defined by 

kg = VgZe (B 2 ) /(cIBg dTi/dlfj) (8) 

so kg — > k\\ in the local limit. That ku ^ kg in the pedestal is a central new result of this 
work. 

As a result of these flow modifications, the current also changes. We write the electron 
distribution f e = f Me exp(e$i/T c ) + h (in the gauge E\\ = B (E\\B) / (B 2 ) - V||<&) with / Me 
the electron Maxwellian. Keeping 0(1) and 0(6) terms in the electron kinetic equation with 
independent variable w = m c v 2 /2 — e$ , (assuming ^/m c /mi <C 8,) 

fMe dh evu(E\\B)B 

v\\V\\h + (v mc + v El ) ■ V/mc + e$i« mc • + ev\\— V H ^ a + '' N " ' / Me = C c . 

l e OWq ±e\tf ) 

(9) 

Here v me is the electron magnetic drift, C c = C cc + C cl is the electron collision operator, 
C ci « v cl L{h} + fM C v c <m c v\\V\\/T c , L = (l/2)(d/d£)(l -£ 2 )(<9/«90, and £ = v n /v. Expanding 
h = h + hi with hi/ha ~ 5, the leading order solution of (i.e. neglecting $! and 
dT[/dip terms) gives ho representing the usual Pfirsch-Schliiter and bootstrap currents but 
without the dTi/dip contribution. At next order, hi = fMefn e v\\ Vi\\/T e — clriihT^dTi/ dip) / e — 
PocI 2 (dn e /dip)(dTi/dip)h$/e where po = Viiriic/ (ZeB av ), B 2 ^ = (B 2 ), and hr { and h$ are the 
solutions of 

Dh Ti = f M em e (n e T c y l {B 2 )' 1 v\\V\\ {v\\Bk\\) , (10) 
aDh® = v E i ■ V/mc + e$iv mc • + ev\\—^-V\\$i 

l e OWq 

with D = U||V|| — C cc — u ci L and a = p d 2 e~ 1 (dn e /dip)(dT i /dip). Applying f d 3 v to (fTDl . 
J d 3 vv\\hT i = olt-B + k\\B / (B 2 ) and J d?vv\\h<$> = a<$>B — n g /(ZB) where aj; and a$ are 
flux functions, n g = Ti(p Irii dTi/dip) -1 J d 3 v g is the 0(1) normalized density perturbation, 



and we have invoked quasineutrality. Then forming j\\ = e J d 3 v(Zfi — f e )v\\, 



. _cl dp ( B 2 \ dn c B dTi 
311 ~ ~Bd^ [jB 2 ) ~ J + ZlB 2 )^ 

Pod 2 dn e dTi ( (n g ) B 



Z dip dip V (B 2 ) 




where p = p c + p- x . The dp/ dip and \j\\B) terms arise in the local case; the former is the 
standard Pfirsch-Schliiter current, and the latter is the Ohmic and bootstrap contribution. 
The fc|| and n g terms however have not been reported previously. Curiously, the n g term is 
quadratic in the gradients. The Ohmic and bootstrap contribution is 

/ . D \ /771 D \ t {£31 d P , £32 dT e £ T . dTi C nT p I dn c dTi\ 

( 3ll B) = a nco (E.B) - cI Pc (— _ + — — - W ^~ ) (12) 



using notation of Ref. |21|, where er nco , £ 31 , and £ 32 are calculated in the standard way, 
and £t; = \B J d 3 v v\\hT^) and C n r = \B J d 3 v v\\h$} are new dimensionless coefficients. 
In the local case of constant fc||, f lTUj) shows Cr i oc k\\. However, to determine £t ; in the 
global case, (jTOj) must be solved accounting for the poloidal variation of k\\. As with the 
flow, the total current is divergence- free: (jlip . ([6]), and quasineutrality imply (after some 
algebra) = V ■ j = V • (j\\b + cB~ 2 B x V , where the ion plus electron stress 

is computed from (J3J) and f c ~ /m c (1 + e$i/T e ). The new k\\ and n g terms in ( ITTi) arise for 
the same reason as the usual Pfirsch-Schliiter current: a parallel return current must flow 
to maintain V • j = given the perpendicular diamagnetic current. In the pedestal, the 
pressure variation on a flux surface becomes sufficient to modify this diamagnetic current. 

We now discuss our numerical method for solving the pedestal ion kinetic equation. The 
radial domain is an annulus containing the pedestal, several pg wide. As r v ,rT ^> pg, we 
take 77 and T{ constant over this domain for simplicity. Also, radial variation of /, B, and 
V\\0 is neglected. We specify rii(ip), which determines $0 = {Ze)~ l Ti \n(r]/n\). On either 
end of the radial domain, ni(ip) and §o(ip) are uniform for several pg, as in figure [TJa-b, 
allowing local solutions to be used for inhomogeneous Dirichlet radial boundary conditions. 
We discretize in the variables (ip, 8, v, £). 

To solve (J4]), dg/dt is first added to the left-hand side, and with the local solution as an 
initial condition, g is evolved to equilibrium using the following operator-splitting method. 



Consider the successive backwards-Euler time steps 

[&+(i/2) - 9t] /At + #nl{&+(i/2)} = 0, (13) 
fe+i - 0t+(i/2)] /At + ^l{^+i} = Ci + 5, (14) 

where i^NL = (^m - Vi/})(d/di/;) v £ is the "nonlocal" term, and in K L = (v^b + v d0 ) ■ (V) Mi w — 
Ci — i^NL) is om y a parameter. In the sum (fT3|) + (fT4l) . g t+ (i/2)/At cancels, leaving an 
equation equivalent to first order in At to a step with the complete operator + K^. 
However, (|T5|) and (I14p are much easier than a step with the total operator because the 
dimensionality is reduced. Also notice the local and nonlocal operators at each grid point 
need only be L?7-factorized once, with the L and U factors reused at each time step for 
rapid implicit solves. 

Our approach to implementing the full Fokker-Planck field operator, similar to the lo- 



221 ] H and G as unknown fields 



cal code in Ref. 0, is to treat the Rosenbluth potentials 
along with g, and to solve a block linear system for three simultaneous equations: (fl4l) . 
VyH = —4:7rg, and Vi^G = 2H , with V|> the velocity-space Laplacian. Our local solver has 
been successfully benchmarked against many analytic formula and against results of another 
Fokker-Planck codeQ]. More details of the numerical implementation will be described in a 
forthcoming publication. 

The heat fluxes at the two radial boundaries are different due to the different densities, so 
heat will accumulate in the simulation domain, precluding equilibrium unless an appropriate 
heat sink is present. In a real pedestal, there will be a divergence of the turbulent fluxes, 
which could act as this sink in the long- wavelength (drift-kinetic) equation we simulate here. 
Determining the phase-space structure of this sink from first principles is beyond the scope 
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of this work, so we use S = —7 + g(—£)) for constant 7, resembling the sink in Ref. 
for global Sf gyrokinetic codes. Varying 7 by several orders of magnitude or using different 
forms of S cause little change to the results. 

Figures [T]{2] show results of the global calculation for a pedestal with e = 0.3, B = 
B /[l + ecos(0)], and V\\8 =constant. The density decreases by 3x from the top of the 
pedestal to the bottom, varying = v^/ (e 3 / 2 v{V\\9) from 1 — 0.3. The electric field profile 
consistent with this density profile for r v ^> pe is shown in figure [Ub. The electric field 
reaches a maximum magnitude of ~ —0.5viBg/c in the middle of the pedestal. In these 
plots, the radial coordinate r/p 9 is defined by r/p e = ZeB (m i cViI)~ 1 ip where B is the 

8 



toroidal field on axis; r = is an arbitrary minor radius, not the magnetic axis. For the 
sink, 7 = 0.1 where uj t = v{V\\6 is the ion transit frequency. The simulation is run 
to t = lOO/wt, since doubling this duration produces negligible difference in the results. 
Figures [TJc-d and |5] show the parallel flow coefficient k\\ and the normalized poloidal flow kg. 
For comparison, the local ku = kg is also shown, computed at each r by numerical solution 
of (J4]) without Vd or S. Even in the local case, fcii and kg vary slightly across the pedestal 
due to the change in collisionality. Outside of the pedestal, as expected, ku and kg computed 
by the global code are equal, constant on each flux surface, and unchanged from the local 
(conventional) result. Inside the pedestal, fen and kg differ from the local result, and both 
coefficients vary poloidally and change sign. The most dramatic change is a well in ku and 
kg at the outboard midplane. Although the distribution for an up-down symmetric B field 
has the symmetry g(—6, —v\\) = —g{6, v\\) in the local case, in the global case the drift terms 
in the kinetic equation break this symmetry, so the global curves in Figure [2] lack definite 9 
parity. To verify mass conservation, the v\\, u, and v m terms in (jHJ) were each independently 
computed from g, and it was verified that the result indeed summed to zero. 

To conclude, in this work we have demonstrated an extension of neoclassical calculations 
to a density pedestal with r n ~ pg but 3> Pe-, retaining effects of finite orbit width, 
collisionality, and aspect ratio. The kinetic equation remains linear, and a 5f approach is 
possible. A numerical scheme was illustrated, demonstrating convergence on a laptop for 
experimentally relevant parameters. The Rosenbluth potentials are solved for along with the 
distribution function at each step, allowing use of the full linearized Fokker-Planck collision 
operator. 

The analytic and numerical calculations show that in a pedestal, the plasma flow can 
differ significantly from the conventional prediction. While the poloidal flow is oc Bg in the 
core, the same is not generally true in the pedestal, and while the numerical coefficients 
in the parallel and poloidal flow are identical in conventional theory, in the pedestal these 
coefficients ku and kg are generally different. These modifications may be important for 



comparisons of experimental pedestal flows to theory. |24j| Two new contributions to mass 
conservation become important which are normally neglected: E x B motion of the per- 
turbed density, and diamagnetic flow of the pressure perturbation. In general, the poloidal 
flow and dTj dip component of the parallel flow can differ in both magnitude and sign relative 
to local theory, as shown in the figures. 




FIG. 1. (Color online) a) Equilibrium density, normalized to its value at the left boundary. As 
happens to be 1 at this boundary and T\ ~ constant over the domain, this plot also gives the 
v* profile, b) Normalized radial electric field —cI(viBo)~ 1 d$>o/d'ilj. c) The dTi/d^-driven parallel 
flow fen computed in the local approximation (dashed curve) differs from the global result (nearly 
indistinguishable solid curves) in the pedestal. The global code is well converged, demonstrated 
by changing each resolution parameter by 2x. d) Normalized poloidal flow kg and /c||, evaluated 
at the outboard (9 = 0) and inboard (9 = ir) midplanes. 



-0.2 
-0.4 







r/p = 0.5, local k = k e 

r/p =0.5, global k 

r/p =0.5, global k 

r/p = -0.25, local k = k 

r/p s = -0.25, global k 

r/p = -0.25, global k 









-7i -n/2 n/2 n 

e 



FIG. 2. (Color online) Poloidal variation of the parallel flow coefficient k\\ and normalized poloidal 
flow kg at two radial locations straddling the pedestal. 



10 



Associated with the modification to the flow, the usual division of the parallel current into 
Pfirsch-Schliiter and Ohmic-bootstrap components is changed (Eq. pip ), and the dTi/dijj 
contribution to the bootstrap current is altered. In the weak-T/ orderings used here, the 
associated terms are necessarily smaller than adjacent dp /dip terms. However, analogous 
modifications to the current would presumably occur in a full-/ calculation when tt ~ p#, 
giving order-unity departures from local theory in that case. 

We are grateful to Peter Catto and Felix Parra for enlightening discussions and for read- 
ing the manuscript. We also thank S. Kai Wong and Vincent Chan for assistance with 
benchmarking our local code to that of Ref. [7|. This work was supported by the Fusion En- 
ergy Postdoctoral Research Program administered by the Oak Ridge Institute for Science 
and Education. 
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